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Abstract 

Background: Hitchhiking mapping and association studies are two popular approaches to map genotypes to 
phenotypes. In this study we combine both approaches to complement their specific strengths and weaknesses, 
resulting in a method with higher statistical power and fewer false positive signals. We applied our approach to 
dairy cattle as they underwent extremely successful selection for milk production traits and since an excellent 
phenotypic record is available. We performed whole genome association tests with a new mixed model approach 
to account for stratification, which we validated via Monte Carlo simulations. Selection signatures were inferred 
with the integrated haplotype score and a locus specific permutation based integrated haplotype score that works 
with a folded frequency spectrum and provides a formal test of signifance to identify selection signatures. 

Results: About 1,600 out of 34,851 SNPs showed signatures of selection and the locus specific permutation based 
integrated haplotype score showed overall good accordance with the whole genome association study. Each 
approach provides distinct information about the genomic regions that influence complex traits. Combining whole 
genome association with hitchhiking mapping yielded two significant loci for the trait protein yield. These regions 
agree well with previous results from other selection signature scans and whole genome association studies in cattle. 

Conclusion: We show that the combination of whole genome association and selection signature mapping based 
on the same SNPs increases the power to detect loci influencing complex traits. The locus specific permutation 
based integrated haplotype score provides a formal test of significance in selection signature mapping. Importantly 
it does not rely on knowledge of ancestral and derived allele states. 
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Background 

Linking genotype to phenotype is one of the central 
questions in biological sciences. Current approaches to 
map intraspecific variation to causative sequence varia- 
tion use either a quantitative genetics framework (asso- 
ciation mapping) or rely on population genetic theory 
(hitchhiking mapping). 

Population genetic theory predicts that a favorably 
selected allele is either lost or increases in frequency 
until fixation [1]. With the spread of a beneficial allele, 
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linked, non-selected sites also increase in frequency, a 
phenomenon that has been termed hitchhiking [1]. 

Based on this principle, genome scans were performed 
in a large number of species such as human, maize, Droso- 
phila, Arabidopsis thaliana and Plasmodium falciparum 
[2-10]. Selection signatures in cattle based on SNP data on 
single chromosomes were reported on Bos taurus (BTA) 
chromosomes 6 [11], 19 [12] and 29 [13]. Barendse et al. 
[13], Gibbs et al. [14] and Hayes et al. [15] published gen- 
ome wide maps of diversifying selection between Bos 
taurus dairy and beef cattle, Flori et al. [16] between three 
different French dairy cattle breeds, and Gautier et al. [17] 
among several West African cattle breeds. Qanbari et al. 
[18] employed an extended haplotype homozygosity test 
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and published a genome wide map of recent selection 
within the German Holstein dairy cattle population. Gau- 
tier et al [19] also used this signature of selection within a 
recently admixed Caribbean cattle breed. Furthermore 
these authors employed a modified version of Rsb scores 
proposed by [20] to detect local excess or deficiency from 
a given ancestry relative to the average genome admixture 
levels. Qanbari et al. [21] recently published a genome 
scan in several dairy and beef breeds including German 
Brown Swiss cattle based on integrated haplotype scores 
and when contrasting breeds employing F ST statistics. 
However, disentangling selection from nuisance signals 
caused by the demographic history of a breed or species 
based on genome wide polymorphism data remains 
challenging. 

Stringent artificial selection resulted in an enormous 
improvement of production traits over the last couple of 
decades, especially for traits with moderate to high herit- 
ability. In combination with the availability of high den- 
sity SNP arrays and high quality phenotypes, this intense 
selection renders the genome of dairy cattle an optimal 
model to look for signatures of recent positive selection. 

While for genetic model organisms very powerful 
genomic tools are available, these species frequently lack 
phenotypic records to link signatures of selection in the 
genome to actual variation in phenotype unless a huge 
additional phenotyping effort is undertaken. This is the 
great advantage of using livestock species, as numerous 
production- and fitness traits are routinely recorded and 
used in breeding value estimation. 

The estimated breeding value (EBV) expresses the 
genetic merit of a breeding animal estimated based on 
their own performance and performances of all available 
relatives. In the case of dairy bulls this typically includes 
hundreds to thousands of daughters. Furthermore EBVs 
are corrected for systematic environmental effects. 
Therefore the breeding value of an animal is the sum of 
its genes' additive effects based on Fishers infinitesimal 
model [22], which assumes a very large (effectively infi- 
nite) number of loci each with very small effect. Although 
only approximatively correct, application of this model in 
selection paved the way for efficient livestock breeding. 

Since Sax's experiments with beans in 1923 [23] we 
know however, that there are so called quantitative trait 
loci (QTL) that have a bigger than infinitesimal effect 
and that these loci can be mapped i.e. via linkage ana- 
lyses. Such QTL mapping studies as a quantitative 
genetics approach have been very successful in cattle, 
see [24-26] for a summary. 

Rapid improvements in high throughput SNP geno- 
typing technologies and commercially available high 
density SNP arrays for livestock species allowed live- 
stock geneticists to turn towards whole genome associa- 
tion (WGA) mapping approaches in the recent past e.g. 



[27-29] or see [30] for a review in livestock. The number 
of individuals that need to be genotyped to achieve rea- 
sonable power in a stand alone WGA is nevertheless 
still limiting [31]. 

Population genetics provides information that is inde- 
pendent of phenotypic information on putative loci 
under strong directional artificial or natural selection. 
We show in this paper that combining a population 
genetics signal with association tests based on quantita- 
tive genetics in a composite statistic, increases power and 
reduces the number of false positive signals for localizing 
the source of selection. 

In a similar vein, [32] proposed a composite test statis- 
tic of several selection signature signals to increase power 
to detect selection. Barendse et al. [13] discussed the 
potential of combining genome wide scans for selection 
and whole genome association studies. However, as these 
authors were looking for signatures of diversifying selec- 
tion based on F ST values the combination with associa- 
tion results is not straightforward. Akey at al. [33] 
followed up a region on dog (Canis familiar is) autosome 
13 that showed evidence for selection in the Shar-Pei dog 
breed with association mapping and finally dissected the 
molecular basis of the typical skin wrinkling phenotype 
in this breed. Ayodo et al. [34] found that in a case- 
control candidate gene approach in humans the statistical 
power to detect disease variants can be increased by 
orders of magnitude by weighting candidates by their evi- 
dence of natural selection. 

Our composite statistic combines a long-range haplo- 
type statistic, based on genomic signatures of (new) 
positive mutations that are not yet fixed in a single 
population, and the regression coefficient based on 
allele-count indicator variables of a WGA - as the quan- 
titative genetic approach. Both estimators rely on the 
underlying linkage disequilibrium (LD) between the cau- 
sal variant and the genotyped SNP. We further propose 
a new mixed model approach to account for stratifica- 
tion in population based association studies, and we 
introduce a modified extended integrated haplotype 
score test statistic to detect selection. Using computer 
simulations and real data we show that the combination 
of both tests increases the power for localizing the target 
of selection relative to a single test and reduces the 
number of false positive signals. 

Methods 

Experimental Design 

The highest selection pressure in the overall breeding 
goal in Brown Swiss cattle over the last decades was put 
on protein yield, the main trait of interest in this study 
to ensure high power for both mapping approaches. 

The 140 highest and 148 lowest bulls with respect to 
protein yield EBV and a minimal EBV-accuracy (r 2 , 



Schwarzenbacher et al. BMC Genomics 2012, 13:48 
http://www.biomedcentral.eom/1 471 -21 64/1 3/48 



Page 3 of 1 9 



degree of determination) of 0.9 were chosen out of 973 
progeny tested Brown Swiss bulls for selective genotyp- 
ing [35]. Up to two generations were present among the 
genotyped bulls. The bulls descend from 90 different 
sires and 121 maternal grandsires. Sire and maternal 
grandsire family size ranged from 1 to 20 and 1 to 34 
members, respectively. 

Phenotypes 

Sire EBVs were obtained from the genetic evaluation 
centre LfL Grub, Germany from the August 2008 
genetic evaluation for PY. EBVs for protein yield are in 
kilogram units. 

Genotypes 

Genomic DNA was prepared from semen straws follow- 
ing standard protocols using proteinase K digestion and 
phenol-chloroform. Across all samples the concentration 
was set to 50 ng/|il. Bulls were genotyped according to 
the manufacturer instructions with the Illumina Bovi- 
neSNP 50K Bead chip® comprising 54,001 SNPs at the 
Institute of Human Genetics of Helmholz Zentrum 
Munchen, Germany. Genotypes of one individual were 
omitted due to a call rate of < 90%. The average call rate 
of the remaining 287 bulls was 98.6% corresponding to 
approximately 53,230 genotypes obtained per individual. 
The software PLINK, version 1.03 [36] was used to filter 
raw genotype data. SNPs with known genomic location 
on autosomes, with a minor allele frequencies of > 5%, 
that were missing in less than 10% of bulls were consid- 
ered. We then filtered for all SNPs for which the ances- 
tral state of the allele was reported by [37]. The final 
dataset contained 34,851 SNPs. Haplotypes were inferred 
with fastPHASE, version 1.2 [38]. Parameters in fas- 
tPHASE were set to 10 random starts for the EM algo- 
rithm and 10 clusters. Haplotypes were inferred for 
whole chromosomes ignoring pedigree information. 

Detection of Selection Signatures 

We wrote R and C++ scripts to calculate extended hap- 
lotype homozygosity (EHH) test statistics from phased 
haplotype data as proposed by [6,7]. Briefly, the EHH of 
a core SNP is calculated as: 



EHH = 



where q is the number of samples of a particular core 
SNP allele i, e^ is the number of samples of a particular 
extended haplotype j, carrying the allele i at the core 
position, and s is the number of unique extended 



haplotypes [7]. EHH captures, as a function of distance, 
the decay of identity of haplotypes that carry a specific 
core allele. EHH starts at one and decays to zero, with 
increasing distance for both alleles at each core SNP. 
The area under the EHH curve that results from plot- 
ting EHH versus distance is expected to be greater for 
the selected allele than for the neutral allele. As pro- 
posed, we computed the integrated EHH (iHH) as the 
integral of the observed EHH in both directions from 
the core position until EHH reaches 0.05 [9]. The 
unstandardized integrated haplotype score (uiHS) is 
then calculated as 



uiHS = In 



/iHH, 



ancestral 



\ ^^^derived 



Voight et al. [9] estimated the expectation and stan- 
dard deviation (SD) of In (iHH ancestra i/iHH derived ) in bins 
of derived allele frequencies from the empirical distribu- 
tion at SNPs whose derived allele frequency p matches 
the frequency at the core SNP. The resulting standar- 
dized iHS (iHS Voight ) follows approximately a standard 
normal distribution. 



iHS Voight 



uiHS — E p (uiHS) 
SD p (uiHS) 



Since standardisation is based on the frequency of the 
derived allele this sets an upper limit to the age of the 
mutation. This test statistic answers the question of how 
unusual the length of a haplotype is, assuming the same 
age of allele across all observed selection coefficients 
acting on any core SNP with a similar derived allele fre- 
quency in the genome. It therefore does not provide a 
formal test of significance. Furthermore if different out- 
groups are used to define ancestral and derived states 
this sets different age boundaries to the mutations 
resulting in less precise standardisation. 

A locus specific permutation-based iHS 

When the rate of EHH decay is similar for the ancestral 
and derived allele, as expected for a neutral locus, uiHS 
is ~ 0 [9]. 

Voight et al. [9] showed via simulation that extremely 
positive and negative iHS scores are both potentially 
interesting selection signals and polarisation with the 
ancestral allele results in a change of sign, but does not 
change the magnitude of the uiHS test statistic. 

In the following we introduce a locus specific permu- 
tation based approach that relies on minor and major 
allele frequencies rather than ancestral and derived 
states, respectively. Most importantly this test statistic 
provides significance of deviations of uiHS from its neu- 
tral expectation. 
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The core site in iHS test statistics is used to define two 
groups of haplotypes for comparison with regard to their 
block structure. We shuffled core SNP alleles 1,000 times 
at the core position while retaining the neighbouring 
haplotype configuration and calculated uiHS for each 
permuted sample (iHS P ) within each core SNP. Random 
shuffling SNP alleles at core sites randomizes allocation 
of haplotypes to the two groups for comparison while 
maintaining the LD structure in the surrounding geno- 
mic region. This simulates the null hypothesis of neutral- 
ity: the site was not subject to selection. We hereby 
obtain an empirical distribution of iHS under the HO for 
each SNP, from which we obtain the probability that we 
see such an extreme iHS just by chance. The locus speci- 
fic standard deviation of the 1,000 iHS P test statistics is 
then used to scale the observed deviation of uiHS from 
its expectation zero. Scaled, permutation based iHS 
(siHSp) is therefore calculated as 

uiHS - 0 

siHSp = 

SD(iHS P ) 

Since the empirical mean of permuted iHS statistics is 
approximately 0 (see Additional file 1, Figure SI) our 
test is a formal test of significance, given the allele fre- 
quency of the core site and the LD structure in the sur- 
rounding region. This is a property of crucial 
importance of a test statistic, especially since we want to 
combine our results with an association test from a 
WGA study. 

Generally SNP sites with low minor allele frequencies 
show larger SD of iHS P . However, to avoid any additional 
bias due to possibly remaining dependence of siHSp on 
allele frequency, siHSp was fit in a linear model by regres- 
sing the SNP minor allele frequency (MAFi) at the core 
site on siHSp. For each site, the random residual was 
obtained and subsequently standardized using the stan- 
dard deviations SDfej) of residuals across all SNPs. In 
contrast to [9] (here termered as iHS Voight ), our fre- 
quency correction is not done based on the expectations 
of SNPs within allele frequency bins but carried out on a 
continuous scale. 

siHSp ij = MAFi + £ij 

The resulting frequency corrected, and scaled test sta- 
tistic is termed iHS. 



This final test statistic is approximately standard nor- 
mally distributed. 

Since no high resolution genetic map was available for 
the SNPs in this study, physical distances between SNPs 
were used for calculating all integrated haplotype scores. 



Whole Genome Association Study 

Standard statistical tests, e.g. regressing phenotype on 
allele count in a linear model, are inappropriate for 
population based WGA in structured populations 
because they either result in an inflated proportion of 
spurious marker - phenotype associations or mask true 
associations (e.g. [39,40]), even with modest levels of 
population stratification and/or admixture. In the case 
of cattle populations, artificial insemination schemes 
allow the use of a few superior bulls as sires of the next 
elite sire generation. Thus, genotyped bulls are fre- 
quently paternal half-sibs or share the same maternal 
grandsire. The subsequent family structure can cause 
substantial stratification. Quantile - quantile plots and 
inflation factors (A) were used to characterize the extent 
to which the observed distribution of P-values follows 
the expected null distribution. Inflation factors were cal- 
culated as 

X = bAedian{T\, T 2 , . . . T^)/0.456 

with T? = PllVar{Pi), where is the effect of the i-th 
SNP (i from 1 to N),Var(/3,) the variance of the estimate 
and 0.456 the median of the Xuf distribution [41]. 

Recently, linear mixed models were proposed to effec- 
tively account for different levels of relatedness by incor- 
porating pairwise genetic relatedness into the model 
[31]. This approach relies on the fact that the pheno- 
types of two genetically related animals are more similar 
than those of genetically distant individuals. Estimation 
of covariance between individuals is assisted by the 
availability of a marker based kinship matrix, which can 
be estimated more accurately using genotype data from 
the WGA experiment than from pedigree information. 

We therefore employed the following single locus 
mixed model which we term "MIX" that explicitly mod- 
els the polygenic relationships among inviduals, as 

y = Xb + Za + e, (1) 

where y is a vector of sire EBVs for protein yield, X is 
the design matrix in which SNP genotypes were coded 
0, 1 and 2, counting the number of minor alleles and b 
the vector of regression coefficients on recoded SNP 
genotypes. Z denotes the design matrix for random 
effects with a ~ N (0, Ga a 2 ) being the vector of poly- 
genic effects, a a 2 the additive genetic variance and G the 
genetic covariance matrix and e ~ N (0, Ia e 2 ), a vector 
of residual effects. G was obtained from pairwise identi- 
cal by descent (IBD) estimates using genome wide SNP 
data as implemented in PLINK [36], in which the IBD 
state is estimated by a hidden Markov model, given the 
observed identity by state (IBS) sharing and genome 
wide levels of relatedness between the pairs. Diagonal 
elements of G were calculated as 1+F, with F being the 
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inbreeding coefficient estimated from SNP data using 
PUNK [36]. 

Mixed models were solved in R (http://www.cran.r-pro- 
ject.org) via direct matrix inversion. Empirical P-values 
were calculated by an adaptive permutation procedure, 
shuffling the vector of genotype codes among pheno- 
types. This does not destroy the relationship between 
IBD status and phenotypes, but breaks up any association 
between SNP genotypes and phenotypes. This leaves LD 
patterns unperturbed and hence does not control for 
stratification. The number of permutations was sequen- 
tially increased up to 1 x 10 6 permutations if the SNP 
indicated association. The empirical P-value was calcu- 
lated as the number of test statistics obtained on per- 
muted sets being greater than or equal to the observed 
test statistic. 

All 34,851 SNPs were tested one after the other for 
association with the protein yield (PY) phenotype. 

As model MIX did not overcome the stratification 
present in our highly structured sample we applied a 
two stage approach. Besides accounting for the relation- 
ship via a mixed model, stratification was accounted for 
by pre-correcting SNP genotype codes for sire and 
maternal grandsire (MGS) differences using the follow- 
ing regression model 



gUj = sirej + MGSj + By, 



(2) 



where gt is the recoded genotype code (0, 1 and 2 for 
1-1, 1-2 and 2-2 allele combinations, respectively with 1 
representing the minor allele), sire is the fixed effect of 
sire i and maternal grandsire the fixed effect of maternal 
grandsire j and Sij~ N (0, Ia e 2 ), the vector of random 
residual effects. Sire- and maternal grandsire families 
smaller than five were merged into one group. 

Residuals e,y were used instead of raw recorded geno- 
types (0, 1 and 2) in the design matrix X of equation 
(1), henceforth termed method "MIXStrat". 

Evaluation of WGA via Monte Carlo Simulation 

The proposed method to account for stratification is 
specific to situations typically observed in intensively 
selected livestock species and populations. We evaluated 
the effectiveness of MIXStrat by Monte Carlo simula- 
tions. Phenotypes, sire- and maternal grandsire family 
structure were taken from the population under consid- 
eration. Genotypes for 287 bulls and 10,000 diallelic 
sites were sampled based on the following procedure: 

First, the allele frequency p of the first allele at a SNP 
was drawn from a uniform distribution, the allele fre- 
quency for the second allele q at this SNP is then given 
by q = 1-p. Two alleles each were sampled for all sires 
and maternal grandsire according to these frequencies. 
Bulls inherited sire and maternal grandsire alleles 



following Mendelian rules. Alleles inherited via the dam 
were sampled corresponding to the population allele fre- 
quencies. This simulates the null model (e.g. no effect of 
the locus on the phenotype) taking into consideration 
the observed population structure. Association was 
tested for, using the models MIX and MIXStrat. 

Rate of False Positives 

The average rate of false positive detections across m = 
100 random repetitions was calculated as 



ap ■■ 



E( pg ir^ 05} )M 



i=i 



10, 000 



Power Analysis 

For true associations the mean genotype values within 
bulls of sire and maternal grandsire are correlated with 
phenotypic family means. This information is not uti- 
lized when genotypes are recoded and will thus reduce 
power. We evaluated the power of MIXStrat relative to 
the power of MIX under the alternative model. This was 
achieved by simulating an additive QTL effect which 
explained 1, 5 and 10% of the EBV variance: 



OiQTL 



'EBV 



QTL, 



SIZE 



2/7(1 -p) 



where a QTL is the allele substitution effect [42], a EB v 
the variance of EBV, QTL SIZE is the size of the effect as 
proportion of a 2 EBV and p the allele frequency of the 
simulated diallelic locus. Power was calculated as 

™/Pe{P<a Bonf }\ 

Power = > — \ m 

^\ 10,000 Y 

i=i x 7 

with a Bon f being the 5% Bonferroni- corrected type I 
error threshold of 2.5 x 10" 5 and m being the number 
of random Monte Carlo repetitions. 

Composite Test Combined Significance Test and False 
Discovery Rate 

We used Stouffer's method [43] to combine P- values 
from the association study with those from the selection 
signature analysis {Pcomb)- 
The test statistic was calculated as 



where Z is the standard normal variable under H 0 , z 
(Pi) is the P - value from test i transformed to Z and k 
is the number of tests that are combined in the test sta- 
tistic. P - Values Pcomb were obtained using the 
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quantile function of the standard normal distribution. 
The tail area based false discovery rate (FDR) was calcu- 
lated from Pcomb values using the R package fdrtool, 
vl.2.5 [44]. Significance was declared if the q value (FDR 
corrected P - Value) was < 0.10. 

Results 

Evaluating the locus specific permutation of the iHS test 
statistic to detect signatures of selection and comparison 
to iHS Voight 

We mapped selection signatures with iHS Voight and our 
newly proposed iHS to detect sites under selection. 

Table 1 shows the expectations and standard devia- 
tions for each of the derived allele frequency bins used 
for the frequency correction to calculate the approxi- 
mately standard normally distributed iHS Voight test sta- 
tistic. Differences in expectations among derived allele 
frequency bins (Table 1) necessitate working with an 
unfolded frequency spectrum for iHS Voight . 

Figure 1 shows that the standard deviation of 1,000 
randomly permuted iHS statistics (iHS P ) is nearly con- 
stant at -0.18 for SNPs with a MAF > 15% but increases 
more than two fold for SNPs with lower MAFs. A simi- 
lar trend can be seen for iHS Voight where SD in the <= 
0.1 and > 0.9 is higher compared to the rest of the 
derived allele frequency bins. 

For SNPs with low minor allele frequencies we found 
a relatively higher proportion of extreme unsealed iHS 
statistics. We postulate that this is due to increased 
rates of false positives, since power simulations by [9] 
and [45] show that iHS Voight is powerful for loci with 
intermediate allele frequencies and that the power of the 
test drops substantially when the selective sweep is close 
to fixation, in other words for SNPs with low MAF. 

Figure 2 shows that the proportions of permutation - 
based iHS signals with a P - Value < 0.001, < 0.005 and 
< 0.01 are relatively smaller for SNPs with low minor 
allele frequency. This tendency cannot be seen as clearly 



Table 1 Means and standard deviations (SD) in defined 
frequency bins for uncorrected integrated haplotype 
score (uiHS) test statistics to calculate iHS Voight . 



Frequency of derived allele 


Mean 


SD 


<= 0.1 


-1.04 


1.04 


0.1 - 0.2 


-0.93 


0.94 


0.2 - 0.3 


-0.73 


0.92 


0.3 - 0.4 


-0.48 


0.92 


0.4 - 0.5 


-0.26 


0.93 


0.5 - 0.6 


-0.06 


0.92 


0.6 - 0.7 


0.16 


0.94 


0.7 - 0.8 


0.39 


0.95 


0.8 - 0.9 


0.65 


0.99 


> 0.9 


0.75 


1.06 



for the iHS Voight test statistic. For SNPs with MAF > 
0.20 we see that our iHS test yields a higher proportion 
of significant loci when compared to traditional 

iHS VoigM 

Figures 3 and 4 show the histograms of iHS Voight and 
iHS, respectively. Figure 5 shows a QQ-plot of iHS and 

iHS Voight 

Figures 3, 4 and 5 all suggest that iHS has increased 
power compared to iHS Voight . 

Our permutation based standardization allows a for- 
mal test against the null hypothesis of neutrality at a 
core SNP (expectation zero). Our standardization is 
against 1000 permuted test statistics at the same locus 
in the same LD background. We therefore do not need 
to define the state of ancestral and derived allele. 

Additional file 1, Figure S2 shows a histogram of 
derived allele frequencies and Additional file 1, Figure 
S3 a histogram of minor allele frequencies of the 34,851 
SNPs used in this study. Additional file 1, Figures S4 
and S5 show histograms of P - Values for iHS Voight and 
iHS, respectively. 

Detection of Selection Signatures in the Brown Swiss 
dairy cattle population 

Manhattan plots for iHS Voight and iHS for each auto- 
some except BTA 6 are shown in Additional file 2, Fig- 
ure S6 - S33 plots A and B. 

Among the 34,851 SNPs tested genome wide 1,710 
and 1,621 SNPs had a test statistics > |1.96| with 
method iHS Voight and iHS, respectively. 

Distribution among chromosomes is remarkably 
uneven: BTA 5, 6, 12, 19 harbor 148, 124, 98, 89 sites, 
respectively which corresponds to 8 - 11% of all investi- 
gated SNP on the corresponding chromosomes that 
show significance applying iHS. On other chromosomes, 
namely BTA 28 and 17 ~ 1% of investigated SNPs exhi- 
bit significant selection signatures. 

The same is true for iHS Voight BTA 5, 6, 12, 16 and 19 
have 171, 131, 148, 136 and 112 SNPs that show an 
iHS Voight test statistic > |1.96| which corresponds to 8 - 
14% of all SNPs on these chromosomes. BTA 7, 25 and 
27 have only around -1% sites with extreme iHS Voight 
test statistics. 

One particularly illustrative example is given by SNP 
Hapmap52798-ss46526455 located in the proximal 
region of BTA 14 at 0.565311 Mb (see Figures 6, 7 and 8 
and Additional file 2, Figure S18). An iHS of 4.13 for this 
SNP with a frequency of 0.2 for allele G exhibits a larger 
area under the EHH curve as compared to allele A and 
was possibly under selection. Interestingly, this SNP is in 
close neighborhood to the well known DGAT1 K232A 
polymorphism, located at 0.444-0.447 Mb, with strong 
effects on milk production traits (e.g. [46-50]). In the 
Brown Swiss (BS) breed the frequency of the K allele is 
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Figure 1 Plot of standard deviations of permuted integrated haplotype scores (iHS) (1 7 000 permutations) versus minor allele 
frequency (MAF). 



very rare with about 2% in the German BS population of the Illumina BovineSNP 50K Bead chip®. Interestingly, 

[48] and fixed for the A allele in the Italian BS population our iHS provided a strong and convincing signal of selec- 

[51]. Most likely allele A was selected for because of its tion, while the iHS Voight (0.81) provides considerably 

milk yield increasing effect while it reduces fat content. weaker support. Hence, this might illustrate the increased 

Note that the DGAT1 K to A mutation itself is not part power of our modified iHS as compared to the iHS Voight . 
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Figure 2 Proportion of SNPs with significant selection signals, relative to all SNPs with a minor allele frequency (MAF) below the 
value given along the x-axis. solid line: permutation based iHS, dashed line: iHS Voight ; symbols: circles, squares and rhombs symbolize SNPs 
with P- values for the corresponding test statistic below 0.001, 0.005 and 0.01, respectively. 
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Figure 5 Quantile - quantile plot of P - Values from selection signature analysis for 34,851 SNPs using our modified iHS (X) and the 
jHS voigth ( # j test statist j CS/ respectively. 



However, there is growing evidence for additional poly- 
morphisms in the DGAT1 gene and its neighborhood 
that cause phenotypic variation for milk production traits 



eg [52,53]. Of particular interest is a QTL mapping study 
in the German- Austrian-Italian BS population [54], that 
reported significant QTL for milk yield and protein 
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Figure 7 Haplotype bifurcation plot of Hapmap52798-ss46526455. The top figure shows the sweeping allele "G" while the bottom figure 
shows allele "A". This figure shows the breakdown of LD from the core SNP with increasing distance in both directions. The core SNP represents 
the root of the diagram. Each SNP represents a node and is an opportunity for further branching. If both alleles of a SNP are present on a 
haplotype the line branches. The thickness of the lines corresponds to the number of samples carrying the haplotype. The length of a branch 
corresponds to the distance between SNPs. 
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percent in the DGAT1 region although all bulls in this 
study were shown to be homozyogous for the p.K232A 
polymorphism [55]. This finding is supported by the 
large SNP effects estimated for fat and protein percent in 
the US - BS population (http://aipl.arsusda.gov/Report_- 
Data/Marker_Effects/marker_effects.cfm?Breed=BS) 
albeit the near fixation of allele A in this breed. So it is 
likely that the selection signal that is picked up by iHS is 
not purely for the DGAT1 p.K232A polymorphism but 
for the proximal region of BTA14 as a whole including 
the VNTR polymorphism in the promoter region of the 
DGAT1 reported by [53]. 

Association Study on PY 

We used 34,851 SNPs that met our stringent quality cri- 
teria and also had the ancestral allele reported in litera- 
ture for association testing. Population stratification was 
accounted for by including IBD estimates from the gen- 
otype data (method MIX). A quantile - quantile plot 
analysis indicated, that this procedure did not suffi- 
ciently account for population stratification in our data- 
set (inflation factor X = 1.34) (Figure 9). 

We therefore developed a new strategy to reduce the 
number of erroneous association signals in our data 
(method MIXStrat). Both the quantile - quantile plot 



(Figure 9) as well as an inflation factor X of 1.02 con- 
firmed that the MIXStrat model successfully controlled 
for spurious results caused by stratification of our sam- 
ple. Nevertheless, we also experienced a drop in power, 
as expected. The SNP with the smallest q value in 
method MIXStrat was 0.5639561 (tail-area based false 
discovery rate (FDR)) calculated with R package fdrtool, 
vl.2.5 [44], corresponding to a nominal P-value of 
4.778973e-04. Note that the flattening out of the P - 
Value curve for method MIX is a consequence of the 
adaptive permutation procedure. 

Evaluation of WGA via Monte Carlo Simulation 

Computer simulations showed that using MIXStrat the 
sample size in this study is sufficient to only detect strong 
effects explaining at least 10% of the phenotypic varia- 
tion. The Monte Carlo simulation did not account for 
LD because conservative significance thresholds using 
Bonferroni correction were used. Nevertheless, it assesses 
the influence of population substructuring in single SNP 
regression whole genome association studies. Our simu- 
lations show clearly that the sire-, paternal grandsire- 
and maternal grandsire structure in dairy cattle popula- 
tions alone can create significant results without any 
association between genotype and phenotype. 
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Figure 9 Quantile - quantile plot from association study on protein yield using model MIX (X) and MIXStrat (O), respectively. 
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Additional file 3, Figure S34 shows a histogram of 
allele substitution effects across all 34,851 SNPs tested. 

Rate of False Positives 

Empirical type I error rate a and the inflation factor A 
using MIXStrat were 0.05 and 0.99, respectively, while 
an a of 0.13 and A 1.35 was observed when method 
MIX was applied. Furthermore the quantile - quantile 
plot in Figure 10 shows clearly that a standard mixed 
model cannot fully account for the stratification present 
in the data, whereas our MIXStrat approach succeeds in 
controlling the type I error rate under the simulated 
null distribution. 

Power Analysis 

As expected, MIXStrat reduced power under the model 
of an existing QTL (Table 2). Power reduction was 
0.069 and 0.10 for QTL that explain 10 and 5% of the 
trait variance, respectively. The relative power loss com- 
pared to MIX analysis of 8.6 (10% QTL size) and 30.6% 
(5% QTL size) indicates, that for large and particularly 
moderate QTL sizes our method leads to substantial 
power reduction. 

As further shown in Table 2 our dataset has sufficient 
statistical power to detect QTL explaining > 10% of the 



variance in EBV. Effect sizes of that magnitude are 
expected to be rare in livestock species [56]. MIXStrat 
without integration of information on selection signa- 
tures has insufficient power to detect loci explaining 
only 1% of the variance. 

Consensus of Selection Signature Signals and Association 
Signals 

A positive iHS value indicates that the minor SNP allele, 
relative to the major allele, is associated with the larger 
integrated EHH statistic and was possibly selected for. 
Likewise the estimated regression coefficient in the associa- 
tion analysis (PMixstrat) represents the estimated increase in 
trait value per additional copy of the minor allele. Thus 
alike signs of iHS test statistics and PMixstrat indicate that 
the SNP is causative by itself or is in LD with a causative 
site that is under positive selection. Opposite signs of iHS 
and PMixstrat ma Y De observed when sites have pleiotropic 
effects and were selected on a different, possibly unob- 
served, trait. Generally one would expect to see a higher 
proportion of like signs as compared to opposite signs and 
a positive correlation coefficient for traits of major eco- 
nomic importance in the selection history of a breed. 

Table 3 shows the correlations between iHS test statis- 
tics and allele substitution effects given by b M ixstrat f° r 
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Figure 10 Quantile - quantile plot for protein yield under the null distribution from method MIX (X) and MIXStrat (O). 
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Table 2 Results from power calculations of the Monte 
Carlo simulation; the underlying models of MIX and 
MIXStrat are described in the "algorithm section" of the 
paper. 



QTL size in EBV variance 


MIX 


power of MIXStrat 


1% 


0.026 


0.005 


5% 


0.347 


0.234 


10% 


0.772 


0.727 



PY based on all 34,851 sites, as a quantitative evaluation 
of accordance. As expected the overall correlations 
among all sites was low. The correlations between allele 
substitution effects and iHS among sites identified to be 
under selection however was substantial with 0.466 
among the top 1% of sites and even higher among the 
top 0.1% of sites. iHS Voight however was uncorrelated 
with top 1% sites and showed a lower correlation of 
0.228 among the top 0.1% sites as compared to iHS. 
This further supports our notion that iHS is an 
improved haplotype based test statistic for identifying 
important loci. 

Combining Signatures of Selection with Association Tests 

Selection signature - and association test statistics were 
moderately correlated (Pearson correlation coefficient was 
0.091 for iHS and -0.005 for iHS Voight ) across all 34,851 
SNPs as the majority of SNPs are not in LD with a causa- 
tive locus and therefore not under selection. This justifies 
treating the two sets of results as independent and using 
Stouffer's method to obtain P- Values {Pcomb) from a 
combined significance test. Figure 11 indicates that combi- 
nation of tests increases power of detection substantially. 
We applied a FDR threshold of 0.10, which corresponds to 
a nominal P-value cut-off of 2.149935e-06. 

Additional file 2, Figures S6 - S33 show Manhattan 
plots for each of the bovine autosomes, combining 
model MIXstrat with iHS Voight (plot C) and MIXstrat 
with iHS (plot D). All Manhattan plots are annotated 
with selection signature signals among the top 5% found 
by [21] applying iHS Voight in windows of 500 kB in BS 
cattle (symbol o) and in any of the other breeds investi- 
gated, symbol (x). All plots are further annotated with 
QTL results reported from whole genome association 



studies in the cattle QTL database "Cattle QTLdb" [57]. 
We downloaded the gff3 file for btau4 at [57,58]http:// 
www.animalgenome.org/cgi-bin/QTLdb/BT/download? 
file=gbpBTAU. QTL positions are annotated at the mid- 
points between start and end position of the reported 
QTL. QTL annotated outside the assembled bovine 
autosomes and in reverse direction (end position further 
distal than start of QTL) were filtered. Capital letters 
summarize QTL trait ontology classes: B for meat (beef) 
traits, E for exterior traits, H for health traits, M for 
milk traits, P for production traits, R for reproduction 
traits as classified at animalgenome.org 

Only QTL annotated from WGA studies were consid- 
ered, because of the large confidence intervals of QTL 
positions from linkage studies. 

Additional file 3, Figure S35 shows a histogram of 
Stouffers P - Values combining whole genome associa- 
tion results with model MIXstrat and iHS Voight while 
Additional file 3, Figure S36 a histogram of Stouffers 
P - Values combining whole genome association results 
with model MIXstrat and iHS. 

Plots A and B in Figure 12 report selection signature 
mapping results applying method iHS Voight and iHS, 
respectively. We see a nice agreement for both test statis- 
tics with the selection signatures reported by [21], in the 
same breed around 60 and 70 Mb. Both plots show an 
additional strong signal for selection in the region 
between 80 - 100 Mb which harbours the well studied 
casein gene family. This becomes evident by the large 
number of annotated WGA results in this region. 
Recently [59] reported a long range haplotype affecting 
protein yield and mastitis susceptibility in Norwegian 
Red cattle that was introgressed from a Swedish Holstein 
bull into Norwegian Red. SNPs in this region almost 
reach significance in the combined approach (plot D), 
which clearly demonstrates the increased power of the 
combined approach, as results from stand alone WGA 
applying model MIXstrat was far from signifance for any 
of the tested SNPs. 

Hayes et al. [11] do not provide a supplemental table 
of iHS test statistics, we could therefore not annotate 
our Manhattan plots with their results. Nevertheless the 
topology of their Manhattan plot for BTA 6 is strikingly 
similar to our results and results reported by [21]. 



Table 3 Pearson correlation coefficients (95% confidence intervals) of different iHS statistics with regression 
coefficients from association study for protein yield. 

Method SNPs with MAF < 10% all SNPs 

all (N= 4387) top 1% |iHS| (N = 42) all (N = 34,851) top 1% |iHS| (N = 349) top 0.1% |iHS| (N = 35) 

iHS 0.045 0.197 0.091 0.466 0.559 

(0.01 6-0.074) (-0.1 05-0.467) (0.080-0.1 01 ) (0.380-0.544) (0.277-0.75 1 ) 

iHS Voight [9] 0.005 0.21 -0.005 0.002 0.228 

(-0.025-0.034) (-0.092-0.31) (-0.0158-0005) (-0.107-0.102) (-0.1 14-0.521) 
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Figure 1 1 Quantile - quantile plot from association study on protein yield using model MIXstrat (O) and combined test of selection 
signature iHS test statistics and whole genome associations with model MIXstrat (X). 



When comparing plot C and D in more detail it 
becomes evident that combining iHS Voight and WGA 
results does not give as good agreement between the 
combined iHS and WGA test. This is supported by the 
lower correlation among the top 1% iHS Voight test statis- 
tics and regression coefficients from WGA (Table 3). 

Discussion 

Mixed model and method to control for stratification 

The pairwise IBD matrix obtained by PLINK [60] based 
on genome wide SNP data most likely underestimated 
the relatedness among bulls because the underlying algo- 
rithm estimates population allele frequencies from a pre- 
sumably unrelated sample. This is supported by the 
observation that the average IBD estimate was exactly 
0.254 between 795 paternal half-sib pairs and not, as 
expected, elevated due to underlying distant relatedness. 
Stich et al. [61] used SPAGEDI software [62] to estimate 
the IBD matrix and noted a similar problem. SPAGEDI 
also assumes that random pairs of individuals are unre- 
lated and assigns them a kinship coefficient of zero. 

The „Q+K" method, proposed by [40], is a mixed model 
with Q, a matrix containing population substructure to 
estimate v, the vector of population effects and the kinship 
matrix K, which allows estimation of polygenic 



background effects based on information on familial relat- 
edness from recent coancestry. The authors claimed 
improved control of the type I and type II error rates over 
other methods. 

Applying method MIX instead of a least squares allelic 
regression substantially reduced the inflation factor X 
from 2.02 to 1.34 for PY. When we extended method 
MIX by Q, the matrix on population substructure based 
on clusters, estimated using the „pairwise population 
concordance" criteria [40], A was further reduced to 
1.16 (data not shown) but still did not control for all of 
the stratification. The here proposed method MIXStrat 
was able to remove stratification (X = 1.02) and proved 
an advantage over method „Q+K". 

The Monte Carlo simulation confirms that the proposed 
MIXStrat approach deals correctly with all stratification in 
the data, as under the simulated HO the observed -log P - 
Value distribution follows their expectation for the dataset 
as highly substructured as dairy cattle. If our two-step 
approach had resulted in an overcorrection we would 
expect to see deflation in the quantile - quantile plot. 

Detection of Selective Sweeps 

Alleles under positive selection increase in frequency in 
a population and leave distinct signatures in the DNA 
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Position (Mb) 

Figure 12 Manhattan plots of chromosome 6. Legend Figure 12: Capital letters denote QTLs reported from whole genome association 
studies (WGA) at [58], summarized as QTL trait ontology classes: B.. meat traits, E... exterior traits, H.. health traits, M.. milk traits, P.. production 
traits, R.. reproduction traits; o annotates a top 5% jHS Voight test statistic as reported in by [21] in windows of 500 kb in Brown Swiss, x in any of 
the other breeds investigated; Plot A: jHS Voight test statistics, blue line: threshold identifying the top 5%; B: iHS test statistics, blue line: threshold 
identifying the top 5%; C: combined iHS Voight and WGA results with model MIXstrat, D: combined iHS and WGA result with model MIXstrat; blue 
line is a at 10% false discovery rate threshold. 
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sequence. One of these population-genetics based signa- 
tures is the increased length of the haplotype carrying 
the advantageous allele [6] which is caused by a rapid 
rise in frequency of the mutated allele. This creates tem- 
porary LD with nearby loci. Extended haplotype homo- 
zygosity statistics [6] contrast this signature between the 
ancestral and the derived allele at each locus. 

The challenge is to determine whether a signature is 
due to selection or to confounding effects of population 
demographic history, such as bottlenecks, population 
expansions and population subdivision or simply due to 
drift in a finite population. Two striking bottlenecks 
were estimated by [63] in data from 14 European and 
African Bos taurus and Bos indicus cattle populations. 
The first and most prominent bottleneck occurred 
roughly 1,500 generations ago, which corresponds well 
with the time of domestication in cattle. The second 
less pronounced bottleneck, which occurred approxi- 
mately 50 - 100 generations ago, is most likely caused 
by breed formation. We therefore expect substantial 
demographic noise in our set of selection signature test 
statistics. Furthermore consequent assortative mating is 
expected to leave signatures in the genome that can 
easily be mistaken as a signature of selection. 

We mapped selection signatures with iHS Voight . Large 
negative values indicate regions in which newly derived 
alleles are increasing in frequency in the population. 
Large positive test statistics advocate so called soft 
sweeps, sweep from standing natural variation where the 
ancestral allele is increasing in frequency for iHS Voight . 
As changes in the selection regime of dairy cattle are 
well documented and make sweeps from standing 
genetic variation likely we believe that it is important to 
consider both extreme positive and negative iHS test 
statistics as potentially interesting regions in the cattle 
genome. We developed a permutation - based extension 
to the iHS statistic proposed by [9] for which there is 
no need to determine the ancestral and derived state of 
the alleles but contrasts minor and major allele. Our 
method obtains locus specific standard deviations of iHS 
in simulating the null hypothesis and contrasting against 
an expectation of zero. Compared to (iHS Voight ) [9] our 
method is more conservative for loci with low minor 
allele frequencies. A higher correlation coefficient 
between our iHS and (3 M ixstrat indicates that this is a 
consequence of a decreased rate of false positive detec- 
tions rather than reduced power. Despite successful 
selection signature scans in cattle we note that protein 
yield is a typical quantitative trait for which selection is 
essentially multigenic and therefore likely to undergo 
simultaneous selective sweeps. Chevin and Hospital [64] 
showed that for quantitative traits selection at specific 
quantitative trait loci may strongly vary in time and 
depend on the genetic background of the trait. This can 



blur the signature of selection and the corresponding 
region will go undetected in a genome scan [64]. Given 
the long generation intervals in cattle the number of 
generations of intense artificial selection is still small 
which could result in weak selection signals for alleles 
with small effects. Selection signature mapping applied 
to livestock with similarly strong selection but shorter 
generation intervals could be even more powerful. 

Method to combine Selection Signatures with Association 
signals 

We propose a novel approach to increase the power to 
detect association signals. In this study the statistical 
power to detect an association signal was quite limited, 
but by combining two independent sources of informa- 
tion for QTL detection in genome wide studies: associa- 
tion and signatures of selection, we were able to 
increase power and to reduce the false positive rate. 
Loci that explain variation in economically important 
traits are likely under selection and will often show 
incomplete selective sweeps. Thus there is a good 
chance to observe extreme iHS values among loci that 
show association. This is supported by the positive cor- 
relation of 0.446 between PMixstrat an d iHS for loci 
among the top 1% iHS test statistics. Although many of 
the associations identified by our method are not yet 
confirmed, the concordance with prior results from 
WGA studies indicates that we were successful in 
detecting interesting loci. Fine mapping of QTL involves 
genotying of many more SNPs in the associated region 
possibly supported by resequencing a subset of extreme 
individuals [65] and is often tedious and costly. Thus it 
is highly desirable to eliminate false positive associations 
prior to further investigations. 

Our combined approach has highest power at inter- 
mediate allele frequencies, as both independent sources 
of information (selection signature mapping and WGA) 
have highest power at intermediate allele frequencies. 
Alleles that are not allowed to go to fixation are either 
likely to be under balancing selection (heterozygote 
advantage) or have pleiotropic effects with positive and 
negative effects for the traits under selection. Such loci 
are not expected to show a signature of recent positive 
selection. WGA, given the same size of effect, will have 
equal power to identify such loci and loci under positive 
selection. 

Conclusion 

The combination of WGA with hitchhiking mapping to 
identify a bona fide set of SNPs for candidate gene stu- 
dies is very promising. We argue that our method 
improves power of QTL detection and reduces type I 
error rate by combining two independent sources of 
information. Our approach can of course be extended to 
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all routinely recorded phenotypes, but for a proof of 
principle we restricted our analyses to PY as this trait 
was under most stringent selection over the last couple 
of decades and the bulls were selectively genotyped for 
PY to increase power for the whole genome association 
study. 

Stratification is a substantial problem in WGA studies, 
particularly when carried out in livestock populations. 
Our MIXStrat approach controls the type I error rate, 
however at the cost of reduced power. 

We accomplished a whole genome hitchhiking map- 
ping study and identified roughly 1,600 SNPs displaying 
selection signatures that show generally good accor- 
dance with effects estimated in the WGA study. Our 
extension to the iHS test statistic proposed by [9] 
resulted in a reduced false positive rate in the MAF 
class < 10%, however, it provides reliable P - Values 
only after extensive Monte Carlo simulations. 

Given the substantial increase in power and the reduc- 
tion in false positive signals we recommend using our 
combined strategy rather than stand alone WGA. This 
is especially important in small populations where it is 
not possible to genotype additional animals. 

Additional material 
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